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Helium Nuclei in Quenched Lattice QCD 
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We present results for the binding energies for He and 3 He nuclei calculated in quenched lattice 
QCD at the lattice spacing of a = 0.128 fm with a heavy quark mass corresponding to = 0.8 
GeV. Enormous computational cost for the nucleus correlation functions is reduced by avoiding 
redundancy of equivalent contractions stemming from permutation symmetry of protons or neutrons 
in the nucleus and various other symmetries. To distinguish a bound state from an attractive 
scattering state, we investigate the volume dependence of the energy difference between the nucleus 
and the free multi-nucleon states by changing the spatial extent of the lattice from 3.1 fm to 12.3 
fm. A finite energy difference left in the infinite spatial volume limit leads to the conclusion that 
the measured ground states are bounded. It is also encouraging that the measured binding energies 
and the experimental ones show the same order of magnitude. 
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The atomic nuclei have been historically treated as col- 
lections of protons and neutrons. The great success of 
the nuclear shell model since 1949 0, Q, explaining the 
nuclear magic numbers and detailed spectroscopy, has 
established that protons and neutrons are very good ef- 
fective degrees of freedom at the nuclear energy scale of 
a few MeV. Nonetheless, 60 years later, we know for cer- 
tain that protons and neutrons are made of quarks and 
gluons whose laws are governed by QCD. It is a great 
challenge to quantitatively understand the structure and 
property of known nuclei based on the first principle of 
QCD. This direct approach will be more important and 
indispensable if we are to extract reliable predictions for 
experimentally unknown nuclei in the neutron rich re- 
gions of the nuclear chart. In this article we address the 
fundamental question in the research in this direction, 
namely the binding energies of nuclei. 

Interacting multi-baryon systems have been investi- 
gated by several studies in lattice QCD. Nucleon-nucleon 
scattering was first studied in quenched QCD 0, Q • This 
work was followed by a partially-quenched mixed action 
simulation in Ref. [5J. Extraction of nuclear force be- 
tween two nucleons has been investigated in quenched 
and 2+1 flavor QCD @, 0, Hi- All these studies assumed 
that the deuteron channel is not bound for the heavy pion 
mass, m^>Q.Z GeV, employed in the calculations. Very 
recently, NPLQCD Collaboration has tried a feasibility 
study of the three-baryon system focusing on the quan- 
tum number of S°S°n. They found the interaction to be 
repulsive Q . So far no evidence supporting bound state 
formation in multi-baryon systems has been observed in 
lattice QCD. In this letter we examine the helium nu- 
clei, He and 3 He, in quenched lattice QCD using a heavy 
quark mass at a single lattice spacing. 

The binding energy AE of the nucleus consisting of 
Npf nucleons with the mass mjv is very tiny compared 



with the mass M of the nucleus: AE/M ~ O(10~ 3 ) with 
AE = Nnuin — M. This causes a complicated situation 
that it is difficult to distinguish the physical binding en- 
ergy from the energy shift due to the finite volume effect 
in the attractive scattering system [l(| . One way to solve 
the problem is to investigate the volume dependence of 
the measured energy shift: In the attractive scattering 
system the energy shift is proportional to 1/L 3 at the 
leading order in the 1/L expansion [n| El, while the 
physical binding energy remains at a finite value at the 
infinite spatial volume limit. In our simulation we choose 
three spatial extents corresponding to 3.1, 6.1 and 12.3 
fm, which are much larger than those employed in current 
numerical simulations so as to provide sufficient room for 
the interacting multi-nucleon system. 

A major computational problem with multi-nucleon 
systems in lattice QCD is a factorially large number of 
Wick contractions of quark-antiquark fields required for 
evaluations of the nucleus correlation functions. A naive 
counting would give (2N P + N n )\(2N„ + N p )\ for a nu- 
cleus composed of N p protons and N n neutrons, which 
quickly becomes prohibitively large beyond thrcc-nuclcon 
system, e.g., 2880 for 3 He and 518400 for He. 

This number, however, contains equivalent contrac- 
tions under the permutation symmetry in terms of the 
protons or the neutrons in the interpolating operator. 
We can reduce the computational cost by avoiding the 
redundancy. In case of the He nucleus which consists of 
the same number of the protons and the neutrons, the 
isospin symmetry also helps us reduce the necessary con- 
tractions. After a scrutiny of the remaining equivalent 
contractions by a computer we find that only 1107 (93) 
contractions are required for the He ( 3 Hc) nucleus corre- 
lation function. We have made a numerical test that the 
result with the reduced contractions reproduces the one 
with the full contractions on a configuration. 
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TABLE I: Number of configurations (iVeonf), number of mea- 
surements on each configuration (N meas ), acceptance rate in 
the HMC algorithm, pion mass (m T ) and micleon mass (mjv). 



L 


N con f 


tncas 


accept. (%) 


[GeV] 


m N [GeV] 


24 


2500 


2 


93 


0.8000(3) 


1.619(2) 


48 


400 


12 


93 


0.7999(4) 


1.617(2) 


96 


200 


12 


68 


0.8002(3) 


1.617(2) 



Another technique to save the computational cost is 
a modular construction of the nucleus correlation func- 
tions. We first make a block of three quark propagators 
where a nucleon operator with zero spatial momentum 
is constructed in the sink time slice. In this procedure 
we can incorporate the permutation symmetry of two up 
(down) quarks in a proton (neutron) sink operator. This 
is a simple trick to calculate 2 Nn contractions simultane- 
ously We also prepare several combinations of the two 
blocks which are useful for the construction of the nucleus 
correlators. 

We carry out calculations on quenched configurations 
generated with the Iwasaki gauge action at (3 = 
2.416 whose lattice spacing is a = 0.128 fm determined 
with ro = 0.49 fm as an input [T3|. We employ the 
HMC algorithm with the Omelyan-Mryglod-Folk integra- 
The step size is chosen to yield reasonable 
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acceptance rate presented in Table [J] We take three lat- 
tice sizes, L 3 xT = 24 3 x 64, 48 3 x 48 and 96 3 x 48, to 
investigate the spatial volume dependence of the energy 
difference between the nucleus and the free multi-nucleon 
states. The physical spatial extents are 3.1, 6.1 and 12.3 
fm, respectively. 

We use the tadpole improved Wilson action with 
csw = 1.378 [3]. Since it becomes harder to obtain 
a reasonable signal-to-noise ratio at lighter quark masses 
for the multi-nucleon system, we employ a heavy quark 
mass at k ~ 0.13482 which gives = 0.8 GeV for the 
pion mass and mjy = 1-6 GeV for the nucleon mass. 
Statistics is increased by repeating the measurement of 
the nucleus correlation functions with the source points 
in different time slices on each configuration. The num- 
bers for the configurations and the measurements on each 
configuration are summarized in Table |TJ We separate 
100 trajectories between each measurement with r = 1 
for the trajectory length. The errors arc estimated by 
the jackknife analysis choosing 200 trajectories for the 
bin size. 

The quark propagators are solved with the periodic 
boundary condition in all the spatial and temporal di- 
rections, and using the exponentially smeared source 
A e -B\ x \ after the Coulomb gauge fixing. On each vol- 
ume we employ two sets of the smearing parameters: 
(A, B) = (0.5, 0.5) and (0.5, 0.1) for L = 24 and (0.5, 0.5) 
and (1.0,0.4) for L = 48 and 96. Effective mass plots 




FIG. 1: Effective mass of He nucleus with Si (circle) and S2 
(square) sources at L = 48 in lattice units. Fit results with 
one standard deviation error band are expressed by solid lines. 



with different sources, which arc shown later, help us 
confirm the ground state of the nucleus. Hereafter the 
first and the second smearing parameter sets are referred 
to as "Si,2 n , respectively. 

The interpolating operator for the proton is defined 
as p a = SabcduafCj^d^u" where C = 7472 and a and 
a, 6, c are the Dirac index and the color indices, respec- 
tively. The neutron operator n a is obtained by replacing 
it" by d" in the proton operator. To save the computa- 
tional cost we use the nonrelativistic quark operator, in 
which the Dirac index is restricted to upper two compo- 
nents. 

The He nucleus has zero total angular momentum 
and positive parity J p = + with the isospin singlet 
1 = 0. We employ the simplest He interpolating oper- 
ator with the zero orbital angular momentum L = 0, 
and hence J — S with S the total spin. Such an 
operator was already given long time ago in Rcf. flil ]. 

He = (xr] - xn) /V2, where \ = ( H 1 — ] + [— H H 

-[+ +] -h + + -] )/2 and x = ( [+- + -] + h + -+] 

+ [+ - -+] + h + +-]-2[+ + ]-2[- - ++] )/\/l2 

with +/— being up/down spin of each nucleon. 77,77 are 
obtained by replacing +/— in \ by p/n for the isospin. 
Each nucleon in the sink operator is projected to have 
zero spatial momentum. We also calculate the correla- 
tion function of the 3 He nucleus whose quantum numbers 
are J p = i + , I = \ and l z = \. We employ the inter- 
polating operator in Ref. (l7| with the zero momentum 
projection on each nucleon in the sink operator. 

Let us first present the He nucleus results. Figure Q] 
shows the effective mass plots of the He nucleus correla- 
tors with the Si t 2 sources on the (6.1 fm) 3 spatial volume. 
We find clear signals up to /j « 12, beyond which statis- 
tical fluctuation dominates. The effective masses with 
the different sources show a reasonable agreement in the 
plateau region. The consistency is also shown in the ex- 
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FIG. 2: Effective energy shift of He nucleus in a convention of 
—AEf 1 with Si (circle) and 52 (square) sources at L = 48 in 
lattice units. Square symbols are slightly shifted to positive 
direction in horizontal axis for clarity. Fit results with one 
standard deviation error band are expressed by solid lines. 



ponential fit results in the plateau region as presented in 
the figure. 

In order to determine the energy shift AEl precisely, 
we define the ratio of the He nucleus correlation function 
divided by the fourth power of the nucleon correlation 
function, R(t) = GHe(t)/(GN(t)) 4 , where Gn c (t) and 
GW(i) are obtained with the same source. The effective 
energy shift is extracted as \n(R(t)/R(t + 1)) = -AEf, 
once the ground states dominate in both correlators. In 
Fig.[2]we present time dependence of —AE L S for the Si y 2 
sources, both of which show negative values beyond the 
error bars in the plateau region of 8 < t < 11. Note that 
this plateau region is reasonably consistent with that for 
the effective mass of the He nucleus correlators in Fig. [TJ 
The signals of — AEf are lost beyond i k 12 because 
of the large fluctuations in the He nucleus correlators. 
We determine AEl by exponential fits of the ratios in 
the plateau region, t = 8 — 12 for S\ and t = 7 — 12 for 
S2 , respectively. We estimate a systematic error of AEl 
from the difference of the central values of the fit results 
with the minimum or maximum time slice changed by 
±1. 

The volume dependence of the energy shift AEl is 
plotted as a function of 1/L 3 in the upper panel of Fig.[3j 
Table |H] summarizes the numerical values of AEl at 
three spatial volumes, where the statistical and system- 
atic errors are presented in the first and second parenthe- 
ses, respectively. The results for the Si^ sources are con- 
sistent within the error bars. We observe little volume de- 
pendence for AEl indicating a bound state, rather than 
the 1/L 3 dependence expected for a scattering state, for 
the ground state in the He channel. 

The physical binding energy AE defined in the infi- 
nite spatial volume limit is extracted by a simultaneous 
fit of the data for the Si,2 sources employing a fit func- 
tion of AE + C/L 3 with AE and C free parameters. The 
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FIG. 3: Spatial volume dependence of — AEl = M — Nntun 
in GeV units for He (upper) and 3 He (lower) nuclei with Si 
(open circle) and S2 (open square) sources. Statistical and 
systematic errors are added in quadrature. Square symbols 
are slightly shifted to positive direction in horizontal axis for 
clarity. Extrapolated results to the infinite spatial volume 
limit (filled circle) and experimental values (star) are also pre- 
sented. 



1/L 3 term is added to allow for contamination of scat- 
tering states. A systematic error is estimated from the 
difference of the central values of the fit results using the 
data with the different fit ranges in the determination of 
AE L - The result for AE is 0.0180(62) in lattice units, 
which is 2.9 a away from zero as shown in Fig. [3l The 
error is evaluated from the statistical and systematic er- 
rors added in quadrature. In the following discussions 
we use the combined error. We also try a pure bound 
state fit allowing for an exponentially small finite size 
correction: AE and AE + de~ C2i with AE and Gi )2 
free parameters. We find all the results are in agreement 
with reasonable values of % 2 . 

Based on these analyses we conclude that the ground 
state of the measured four-nucleon system is bounded. 
An encouraging finding is that AE = 27.7(9.6) MeV 
with a -1 = 1.54 GeV agrees with the experimental value 
of 28.3 MeV. However, we do not intend to stress the 
consistency because our calculation is performed at the 
unphysically heavy pion mass, = 0.8 GeV, and the 
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TABLE II: Binding energies of He and 3 He nuclei on each 
spatial volume. Extrapolated results to the infinite spatial 
volume limit are also presented. The first and second errors 
are statistical and systematic, respectively. 



L 




AE L 


[MeV] 






He(Si) 


He(5* 2 ) 


3 He(Si) 


3 He(S* 2 ) 


24 


28(14)(11) 


46.8(7.3)(1.6) 


19.0(6.3)(6.0) 


23.2(3.2)(0.5) 


48 


27(14)(05) 


36(12)(04) 


16.6(6.9)(3.2) 


19.5(5.6)(2.3) 


96 


24(18)(12) 


24(14) (03) 


19.0(7.6)(4.9) 


18.4(6.1)(1.9) 


oo 


27.7(7.8)(5.5) 


18.2(3.5)(2.9) 



electromagnetic interactions and the isospin symmetry 
breaking effects are neglected. 

We also calculate AEl for the 3 He nucleus with the 
5i,2 sources, whose results are presented in Fig. [3] and 
Table ILT1 The trend of the volume dependence is similar 
to the He nucleus case. A simultaneous fit of the data for 
the Si t 2 sources with a fit function of AE + C/L 3 yields 
a finite value of AE = 18.2(4.5) MeV, which means the 
existence of a bound state in the 3 He nucleus channel. 
Our result for AE is about twice larger than the exper- 
imental value of 7.72 MeV. A main reason could be the 
heavy pion mass employed in this calculation. 

As an alternative way to view this result, we com- 
pare the binding energies normalized by the atomic num- 
ber: AE/N N = 6.9(2.4) MeV and 6.1(1.5) MeV for the 
He and 3 He nuclei, respectively. At our unphysically 
heavy pion mass, the three and four nucleon system does 
not show the experimental feature that the binding is 
stronger for He than for 3 He. 

We have addressed the issue of nuclear binding for the 
He and 3 He nuclei. We have shown that the current com- 
putational techniques and resources allow us to tackle 
this issue. Albeit in quenched QCD and for unphysically 
heavy pion mass, we are able to extract evidence for the 
bound state nature of the ground state and the binding 
energies for these nuclei. 

A future direction of primary importance is to inves- 
tigate the quark mass dependence of the binding ener- 
gies of the nuclei. There are several model studies of 
the quark mass dependence of the nuclear binding en- 
ergies [TU which suggest that the quark masses play an 
essential role in a quantitative understanding of the bind- 
ing energies. Another important issue is development of a 
strategy to calculate nuclei with larger atomic numbers. 



The required number of the Wick contractions quickly 
diverges as the atomic number increases, even if the re- 
dundancies are removed with various symmetries. We 
leave it to future work. 

Numerical calculations for the present work have been 
carried out on the HA8000 cluster system at Informa- 
tion Technology Center of the University of Tokyo and 
on the PACS-CS computer under the "Interdisciplinary 
Computational Science Program" of Center for Compu- 
tational Sciences, University of Tsukuba. We thank our 
colleagues in the PACS-CS Collaboration for helpful dis- 
cussions and providing us the code used in this work. 
This work is supported in part by Grants-in-Aid for Sci- 
entific Research from the Ministry of Education, Cul- 
ture, Sports, Science and Technology (Nos. 18104005, 
18540250, 20105002, 21105501). 
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